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Abstract 

Williams and Bjerknes proposed a simple stochastic growth model to describe the 
tumor growth in the basal layer of an epithelium. In this work we generalize this 
model by including the possibility of saturation in the tumor growth as it is clinically 
observed. The time evolution of the average number of tumor cells and its variance 
for both the original and extended models are studied by analytical methods and 
numerical simulations. The generated growth patterns can be compact, connected 
or disconnected depending on the model parameters used, and their geometrical 
properties are characterized through the gyration radius, the number of interfacial 
cells and the density of empty sites inside the patterns. 
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1 Introduction 



The population dynamics, including the growth of normal and tumor cells, 
is a traditional problem investigated in Physics and Biology [1]. In their orig- 
inal work, Williams and Bjerknes [2] built a model (WB model) to describe 
the tumor growth in the basal layer of an epithelium. In their model, one 
phase (the tumor cells) grows faster than the other (normal cells) by a factor 
K, which represents the carcinogenic advantage. This model exhibits two dis- 
tinct behaviors: unrestricted growth (k > 1) and complete regression of the 
tumor {k < 1). In the special 1 the tumor always disappears due 

to the fluctuations. It is worthwhile to mention that the WB model can also 

* Corresponding author. 

Email address: silviojr@fisica.ufmg.br (S. C. Ferreira Junior). 



Preprint submitted to Elsevier Science 



1 February 2008 



be applied to describe the growth of many systems involving two competing 
phases. Using an adequate time step, the WB model can be interpreted as a 
' birth and death' process [3] with constant division and death rates. Thus, all 
the well known results for stochastic processes can be used to understand this 
model. 

However, real tumors exhibit, in addition to the two distinct behaviors pre- 
viously described, a quiescent state in which the tumor size remains constant 
for a long time [4]. The underlying aspects of tumoral biodynamic diversity 
involve a complex set of biochemical processes and environmental constraints 
such as local nutrient availability, mechanical stress, immune response, etc. [4]. 
Numerous models of cancer growth have been recently proposed in order to 
describe the tumor progression, providing auxiliary tools for cancer diagnosis 
and therapy [5]. 

In this paper, we are proposing a generalization of the WB model in which 
the division and death rates of tumor cells depend on their total number. 
Specifically, the probability of cell division decreases whereas the cell death 
probability increases as the number of cancer cells increases. The paper is or- 
ganized as follows. In section 2, the original and the extended WB models are 
defined. Sections 3 and 4 are dedicated to the analytical study of both discrete 
and continuous time versions of these models. In section 5, the geometrical 
properties of the growth patterns generated by the extended model are char- 
acterized through computer simulations. Finally, some conclusions are drawn 
in section 6. 



2 The extended Williams and Bjerknes models 



In the original WB model [2] the tissue is represented by a two-dimensional 
lattice where occupied sites represent tumor and empty sites normal cells. All 
the sites are initially empty, except the center of the lattice, since the tumor 
grows from a single malignant cell, in agreement with the theory of clonal 
origin of cancer [6]. The interfacial cells are defined as those that have one 
or more nearest-neighbor sites of the opposite type. The growth rule is very 
simple: one of the bonds between two opposite cell types is chosen at random 
with equal probability; the normal cell of this chosen bond is replaced by a 
tumor cell, with probability g, or the tumor cell is replaced by a normal one, 
with the complementary probability r — 1 — g. In terms of the carcinogenic 
advantage k, these probabilities can be written as: 



2 



and 



K+1 



(2) 



The limit k = oo corresponds to the Eden model [7], more specifically, the 
Eden model type B according to the definitions given by Vicsek [8]. 

There are many variations of the WB model [9] , but we shall consider only 
the cases that exclude the steps that do not change the pattern configuration. 
Therefore, in all time steps the number of tumor cells increases or decreases 
by an unity with probabilities g or r, respectively. Particularly, we shall study 
the foUwing variation of the original WB model: at each time step a cell 
type is chosen, either a tumor cell, with probability r, or a normal cell, with 
probability g. Then, the selected cell is converted to its opposite type. In other 
words, at each time step, either the division or death of a single tumor cell 
occurs, with probabilities g oi r, respectively. 

Now, we introduce a new model by assuming that the division and death 
probabilities depend on the total number of tumor cells n through Michaelis- 
Mentem functions [1] : 

g(n) = 1 - 3) 

^ r + n ^ ^ 



and 



r{n) — 



an 



r + n 



(4) 



Here < a < 1 and F > are parameters that control the shape of the 
curves. These functional forms were used because they are the simplest ones 
which varies monotonically with n and satisfy the normalization condition 
g{n) + r(n) = 1. These curves are illustrated in figure 1. 



3 Analysis of WB models using stochastic methods 



The stochastic growth rules used in the WB models involve probabilities 
g and r that are explicitly time independent. In consequence, as is the case 
for Markov chains, the chance of these models generate a given pattern in a 
certain time depends only on the probabilities associated to the configurations 
at the previous time. So, the WB growth processes will be described within 
the probability transition equation or master equation frameworks, according 
to the discrete or continuous character of the time [3,10]. 
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Fig. 1. Division and death probabilities in the extend WB model. 
The transition probability equation is written as: 

P{n,t + l)^Y.Tn,mPim,t), 



(5) 



where T^_„ are the elements of the transition matrix. T^^„ provides the transi- 
tion probability from a state containing m tumor cells to a state with n tumor 
cells in the next time. A continuous version of equation (5) is given by the 
master equation: 



^P{n, t)=Yl {Wn,mP{m, t) - Wm,nP{n, t)} , 



(6) 



where Wn,m is interpreted as a probability by unit time or transition rate. 

Initially, the discrete time WB model will be considered. At every random 
selection of an interfacial site to implement an action, the time is incremented 
by an unity. In this way, the transition matrix is: 



T = < 



g ii n — m + 1 
r \i n — m — 1 
if In — ml > 1 



(7) 



This expression is valid only for n > 2, since the WB model is a special type 
of one step process [3] with an absorbent state at n = 0. Indeed, To,i = r and 
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Tifi = 0. Substituting the transition matrix (7) in the probabihty transition 
equation (5) we have: 

P{n, t + 1) = gP{n -l,t)+ rP{n + 1, t) if n > 2 

(8) 

P{n,t+1) =rP{n+l,t) if n = 0, 1 



In this section, we are interested in quantities that depend only on the 
number of tumor cells and not on the spatial distribution of theses cells on 
the tissue. So, the time dependence of the average number of tumor cells {n{t)) 
and its variance cr(t), defined as: 

oo 

{n{t)) = J2nP{n,t) (9) 

n=0 

and 

a'{t) ^ (n'it)) - (10) 
are calculated. Clearly, 

oo 

{n\t))^Y.^'PM. (11) 

n=0 



Using equation (8) and iterating the expressions for < n{t) > and < n^(t) > 
we find: 

{n{t)) ^ n{0) + (12) 



and 



1 - 



K 



K+1 



t. 



(13) 



Equations (12) and (13) show that {n(t)) decreases linearly with time if k, < 
1, and increases linearly if k > 1. In turn, the variance increases with the 
square root of time. Thus, for all k > 1 there is a non-vanishing probability 
of unlimited tumor growth. But, if k = 1, {n{t)) is constant, the variance 
increases with the square root of time and, therefore, independently of the 
initial population, the absorbent state n = will be reached. Moreover, we 
have the maximum variance exactly at k = 1. These exact results found for 
the discrete WB model are also valid for the continuous time approach. 
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4 The extended model 



In this section, the previously described generalized WB model (equations 
(3) and (4)), which includes the possibihty of growth hmitation, is studied by 
a continuous approach based on the master equation [3,10]. Again, the time 
step is defined as the birth or death of a single tumor cell. 



4-1 The macroscopic equation 



First of all, we define the transfer matrix Wn 



Wr,. 



g{m) if n = m + 1 
r{m) if n = m — 1 
if In — ml > 1 



(14) 



or 



Wn,m — g{m)d{m — n — 1) + r{m)5{m — n + 1). 



(15) 



The master equation is obtained substituting the expression (15) in (6): 
d 



dt 



P(n, t) = g{n - l)P{n -l,t)+ r{n + 1) x P{n + l,t)- P{n, t). (16) 



Using (9), (11) and (16) we have: 



d{n{t)) 
dt 



1 - 2a 



n{t) 



T + n{t) 



(17) 



and 



d{n\t)) 
dt 



= 1 + 2 



' n{t)[T + {1 - 2a)n{t)]' 
. r + n{t) , 



(18) 



Notice that (17) and (18) are non-linear equations which cannot be solved 
by analytical or numerical methods. Using a mean field approximation, the 
macroscopic equation for (17) is obtained replacing n{t) by a deterministic 
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Fig. 2. Mean field phase diagram for the asymptotical value of the average number of 
tumor cells. Three possible regimes are found: unrestricted growth, limited growth 
and total regression of the tumor. 

function A^(t). The corresponding macroscopic equations for the first and sec- 
ond moments are: 

^^'(^) , o iV(t)[r + (l-2a)iV(t)] 

dt " + r + iv(t) • ^^^^ 



The star is used in order to distinguish between iV^ and iV^, the second mo- 
ment. 

As one can see, equation (19) exhibits three distinct asymptotic behaviors: 

• A^(t) increases without limit if C(t) > V t ; 

• N\t) ^ if C(t) < V t; 

• A^(t) goes to a non-vanishing constant value if C(t) = at a certain time. 

The phase diagram in the parameter space (F, a) is shown in figure 2. 
The third behavior requires a stable solution [3], which always exists since 
the macroscopic equation is a first order differential equation. This analysis 
reveals the existence of a well-defined phase transition in a = 1/2, i. e., for 
a < 1/2 the average number of tumor cells grows without limit, whereas for 
a > 1/2 this number reaches a constant value i7 > 0. 

In the region of unrestricted growth [a < 1/2), N{t) grows linearly and 
the variance increases as the square root of the time, in agreement with (19) 
and (20) for very large N. Therefore, the model behaves asymptotically as 
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Fig. 3. Numerical solutions of the macroscopic equation for the extended WB model 
around the critical point a = 1/2 and fixed parameter F = 20. Also, the saturation 
value and the crossover time evaluated for these parameters are indicated. 

the original WB model. The transition line a = 1/2 must be considered in 
separated. Fortunately, the equation (19) has a closed solution for a = 1/2, 
Ni{t), given by: 

m{t) = ^{n{o) + ry + 2rt-r. (21) 



Clearly, A^i has the asymptotical behavior Ni{t) = ^/2Tt. Substituting (21) 
in (20), the asymptotical macroscopic approximation to A^^ is obtained and, 
using this result, we find a — y/i. However, numerical simulations suggest that 
(7 = ^ft/^, and this difference will be explained in the subsection 4.3 by taking 
into account corrections to the macroscopic equation. 

In the region of limited growth, the solution Ns(t) reaches a saturation value 
fl for long times. This value can be find by taking ({t) = 0: 

Q = ^ . (22) 
2a - 1 ^ ^ 



Also, the saturation crossover time, (tx), is evaluated through an expansion at 
long times oi Ns{t) around fi. Substituting Ns{t) = il+n{t), where |Ai(t)| -C 
and expanding (19) we have: 



d/i (2a - 1) 



2 



dt 2ar 



pi + 0{pi'). (23) 
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At first order in /i the solution is an exponential decay ii{t) ~ exp{—t/tx)- 
Thus, the saturation crossover time is given by: 



In figure 3 the numerical integrations of (19) around the critical parameter 
a — 1/2 and fixed F = 20 are shown. 

The phase diagram exhibits two lines in the complete regression region. The 
line r = 2a — 1 is the border of the region where N{t) — ^ 0, according to the 
solution of the macroscopic equation. However, if the fiuctuations are larger 
than the saturation value Q, all tumor cells die. Thus, the phase diagram 
shown in figure 2 already take into account the correction considering the 
variance, as evaluated in next subsection. The result is that the region of total 
regression correspond to T < a. 

In order to test the predictions of the macroscopic equation approxima- 
tion, comparisons with simulational results were done. It was observed a good 
agreement between analytical and simulational results for a < 1/2. However, 
for a > 1/2, small differences, that become meaningful when the value of F is 
not large, emerge. All these comparisons are shown in figure 4. As one can see, 
there is a fixed difference between the numerical and simulated solutions in 
the limited growth region. This difference can be found making an expansion 
of the master equation [3] , also called fl expansion. 



4-2 The n expansion 



The Q expansion consists in expand the master equation in powers of the 
characteristic size of the system. In our case, this characteristic size is the Q in 
(22), since it represents an upper bound to the size of the tumor population. 
Indeed, this was the reason for denote the saturation value as Q. 

We assume that the probabilities P„ have a sharp maximum around the 
macroscopic solution with a width of order O^/^. This hypothesis is explicitly 
used in the change of variable: 

n = Q0(i) + Q^/^C, (25) 



where ^ is a new random variable and (f)(t) is a deterministic function. Also, 
the expansion assumes that the rates W{n\m) = Wn,m can be written in the 
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Fig. 4. Comparison between analytical and simulational results. The continuous 
and smooth curves are obtained from numerical integrations of the macroscopic 
equation whereas the other ones axe given by simulations. 2000 samples were used 
in the simulations with the parameters (a) F = 20 and a = 0.49, (b)F = 20 and 
a = 1/2 and (c) F = 50 and a = 0.60. 



form: 



where / and <l>j are arbitrary functions and s = n — m. The technical details 
involved in the expansion can be found in [3]. 

Introducing the notation used in [3] : 
'PuA^) = / s'^xix; s)ds, (27) 
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the expansion up to orders of Q^/^ and DP results in the following equations 
for 0, iO and ): 

^ - ViM, (28) 

^ = <^i,o(0)(e) (29) 

and 

^-M,om^") + ^^M, (30) 



with T = /(Q)/Q. (28)is the macroscopic equation for the system divided by 
fl, (29) gives the correction for (n(t)) in order of fi^/^, and (30) represents the 
first approximation for the fluctuations around the mean. The initial condi- 
tions for (29) and (30) arc (^(0)) = (C^(0)) = 0, because P„(0) = S{n - no), 
i.e., at t = the population is Uq. 

For the extended model, the transition rates (15) can be rewritten as: 



where p = m/Q. In agreement with the notation of (26), we have f{D) = 1, 
$j = if i 7^ 0, and $0 is defined as the right hand side of (31). 

Here, our goal is calculate the asymptotical corrections in time and the 
differences shown in figure 4. Solving (29) and (30) using (31), we find that: 

(e(r))~exp(-^r)^-^0 (32) 



and 



«^W>-55^- (33) 



Therefore, the correction for (n(oo)) given by (32) vanishes, because n = 
il(p{t) +Q^^'^^{t). However, (33) shows that the first correction for the variance 
reaches a constant value: 

-'-»(eirr-J^,. (34) 
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In figure 4(c) the variance was measured as (Tmeasured — 27.5, whereas the ana- 
lytical value calculated by (34) is (Tanaiyticai — 27.4. Thus, there is an excellent 
agreement between the analytical and simulational results for the variance. 
Since the fluctuations reach a constant value which is smaller than the satura- 
tion one, the tumor certainly does not disappear as in the original WB model 
with K — 1. 

The difference between the saturation values obtained by the simulations 
and calculated through the macroscopic equation can not be explained taking 
into account only terms of order Vt^. The first non- vanishing correction in (^) 
involves terms of order 0~^/^[3]: 

= w^y '''' 



Consequently, the first correction in (n) is: 

A = — ^ -. (36) 

2(2a - 1) ^ ' 

Using (36), the calculated correction is ^analytical = 2.5 and the measured 
value in figure 4(c) is /^measured — 2.51, which are in excellent agreement. 
Several other values for the model parameters were tested and a very good 
agreement between simulational and analytical results was found. 



4-3 Asymptotical corrections for a = 1/2 



The expansion cannot be applied for systems without a characteristic 
size. This is the case for equation (21) which, for a = 1/2, gives a solution 
that grows without limit and, consequently, one can not define the character- 
istic size of the tumor. However, it is possible to do an expansion around the 
macroscopic solution (21). Taking n{t) = N^{t)+e{t), where e is the new ran- 
dom variable satisfying \e{t)\/Ni(t) ^ 1, and conserving terms up to second 
order of s/Ni, we find the following equation for s{t): 

d<^> _ _rl^ + + O(e') (37) 



where g{t) = y n(0) + 2Tt. In order to solve (37), it is necessary to know the 
relationship between the first and the second moments of £. Equation (20) gives 
N^{t) — N\ {t)-\-t and the xpansion of — {Ni+ey with the approximation < 
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r? >Ri A^^ results in < >^ t — 2N\ < e >. Substituting this approximated 
relation for < £^ > in equation (37), and taking the asymptotical limit for 
we obtain: 

d < £ > 3 1 
Its exact solution is: 



Considering the correction given by (39), the analytical and simulational 

curves in figure 4(b) become indistinguishable. Now, one can use the result 
(39) to find the correction to the variance. It is possible to show that the 
asymptotical value of the variance, consistent with the approximations used 
above, is given by: 



Equation (40) is able to explain the factor 1/2 present at the simulations shown 
in figure 4(b). Again, the simulational and analytical results are in very good 
agreement for a large number of parameter sets. 



5 Geometrical patterns 



The previous sections were dedicated to the analytical study of the origi- 
nal WB model and its extended version. Now, in this section the simulational 
results focusing the geometrical properties of the patterns generated by the 
extended WB model are presented. The patterns associated to the original 
WB model are spherical and compact with a rough surface. For the extended 
model, the patterns exhibit three distinct morphologies: compact with a rough 
border, connected with internal holes and disconnected with cells isolated from 
each other. These morphologies are shown in figure 5. In the region of unlim- 
ited growth, i. e., a < 1/2, the patterns become compact. This compaction 
occurs because the division probability is always higher than the death prob- 
ability and, as a consequence, all the internal empty sites will be occupied at 
sufficiently long time. The disconnected patterns appear in the region of lim- 
ited growth (a > 1/2). The reason is that the division and death probabilities 
reach the same value p = 1/2, and the cells behave like non-directed random 
walkers. Finally, the connected patterns are generated just on the transition 
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Fig. 5. Typical patterns generated by the extended WB model: (a) compact, (b) 
connected and (c) disconnected. For these patterns F = 1000 were used, a = 0.45, 
0.50 and 0.51 was used in (a), (b) and (c), respectively. The simulations were done 
in lattices with 1200 x 1200 sites and stopped if a tumor cell reaches the border of 
the lattice. 

line a = 1/2. The difference between connected and disconnected patterns is 
that the former does not has a percolation cluster of empty sites inside the 
pattern, whereas in the last this cluster is found. 

These patterns were characterized by its gyration radius (Rg) and the num- 
ber of interface tumor cells (S). The gyration radius is defined by: 



where the sum extends over all the tumor cells and fcm is the mass center 
of the pattern. In the original WB model, Rg and S scale asymptotically 



n 




(41) 



14 



Table 1 

Summary of the exponents found for the extended WB model. The exponents f, a, 
and 7 are defined hy Rg ^ n'^ ,S ^ n"' and p ^ f^. 





a < 1/2 


a = 1/2 


a = 0.51 


r 


I' 


a 




a 


7 


I' 


a 


7 


10^ 


0.50 


0.50 


0.64 


1 


0.18 


0.62 


1 


0.64 


103 


0.50 


0.50 


0.60 


1 


0.18 


0.60 


1 


0.64 


10^ 


0.50 


0.50 


0.54 


1 


0.25 


0.52 


1 


0.63 



with the square root of the number of tumor cells [11]. In turn, only the 
compact patterns {a < 0.5) of the extended model show the same asymptotical 
behaviour for Rg and S. For a > 0.5 the scaling laws change and are dependent 
on the r parameter. Indeed, S* ~ n and Rg ~ n'^, with /^(F) > 1/2, indicating 
that the patterns are fractals with dimensions df = l/u [8]. Notice that, in the 
limited growth regime, n is bounded and the power law is defined before the 
growth saturation. However, since the cells become progessively more distant 
from each other, Rg grows continuously with the time as a power law. The 
numerical results are summarized in table 1. 

The transition from compact to disconnected patterns was studied through 
the density p of internal holes in the patterns, p was defined as the ratio 
between the number of empty and occupied sites enclosed by a circle of radius 
Rg. This definition was used in order to discard the tumor border, where the 
compaction never occurs. For a < 1/2, p decreases, vanishing at sufficient long 
time. In contrast, for a > 1/2, p increases asymptocally as a power law p ~ P. 

In order to test if the exponent 7 varies with the parameter F, simulations 
were done using three distincts F values (10^,10^,10^). For F = 10^ the density 
of internal holes is lower than that for F = 10^, but the corresponding power 
law exponent is higher. Since for both F = 10^ and F = 10^ the simulations 
provide essentially the same values for 7, we suppose that this exponent is 
independent of the parameter F. Thus, we believe that the curve obtained for 
F = 10^ might be a transient behavior. 

The patterns observed in the extended model are very similar to the mor- 
phologies generated by a growth model for primary cancer recently proposed 
[12]. This model considers a complex interaction network among tumor cells 
mediated by growth factors and, in contrast to the present model, involves a 
large number of parameters. Our results are more robust than those observed 
in [12], since the two parameters used in the extended WB model are not 
associated to the properties of the tumor microenvironment. Indeed, they are 
related to macroscopic quantities such as the saturation size of the tumor. 
However, the extended WB model has a limitation. Every compact pattern 
has a unlimited growth and the disconnected patterns always cease their pop- 
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ulational growth. The reason is that the growth only occurs on the interface 
between tumor and normal tissues, whereas in real tumors the internal cells 
also divide. This limitation does not compromise the analytical results since 
the spatial distribution of the tumor cells in the tissue is not considered. 



6 Conclusions 

In this work we propose an extended version of the Williams and Bjerknes 
(WB) model initially used to describe the tumor growth. In this extended 
model, the division and death probabilities correspond to monotonically in- 
creasing and decreasing functions of total number of tumor cells n, respec- 
tively. The average values of n for the original and the extended WB models 
as well as their variances were exactly calculated using stochastic processes 
methods. The original model reveals two possibilities for < n >: unlimited 
growth or total regression of the tumor. However, real tumors exhibit three 
possible behaviors which arc observed also in the extended WB model: the 
ones previously described and a quiescent state where the tumor size remains 
constant for a long time. The differences observed between simulated and 
mean-field results for the extended model were calculated using expansions of 
the involved equations. In addition, the geometry of the growth patterns gen- 
erated by the extended WB model were analysed. Three distinct morphologies 
have been observed: compact, connected and disconnected. All the patterns 
were characterized by its gyration radius, Rg, and number of tumor cells on 
the interface, S. For the compact patterns, Rg ~ y/n and S ~ y/n. For the 
connected and disconnected patterns Rg ~ n" and S ^ n, with v > 1/2, 
indicating that these patterns are fractals. 
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